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Abstract 

This is the second in a series of articles on the numerical solution of 
Friedrich's conformal field equations for Einstein's theory of gravity. We 
will discuss in this paper the numerical methods used to solve the system 
of evolution equations obtained from the conformal field equations. In 
particular we discuss in detail the choice of gauge source functions and 
the treatment of the boundaries. Of particular importance is the process 
of "radiation extraction" which can be performed in a straightforward way 
in the present formalism. 

1 Introduction 

In article [1] we have presented the conformal field equations explicitly in a 
form suitable for solving them numerically. In a brief summary, we have de- 
rived the conformal field equations in the space spinor formalism which is well 
suited to perform the 3 + 1 split because the evolution equations come out in 
a symmetric hyperbolic form (under appropriate assumptions on the free gauge 
source functions) almost automatically. We have also discussed in [1] the fur- 
ther assumption of a hypcrsurface orthogonal symmetry which has been made 
to simplify the implementation. Fixed points of a continuous symmetry usually 
lead to coordinate singularities which have to be treated specially in any finite 
difference method. Therefore, we followed a suggestion of B. G. Schmidt [2] to 
require that there be no fixed points which has the unphysical consequence that 
the global topology of space-time is x R^. However, since the emphasis of 
this project lies in studying the effectiveness of radiation extraction from the 
numerically generated space-times and since these are local methods this is not 
a serious disadvantage. 

In section 2 we first present the numerical method, a Lax-Wendroff method 
in two dimensions and the procedure for choosing the time-step dynamically in 
order to enforce the CFL condition for stability of the algorithm. In section 
3 we show how the boundary can be treated. This is an essential part of any 
numerical scheme because if the boundary conditions are non-physical one has 
to live with the fact that the numerical solution probably differs quite signif- 
icantly in the domain of influence of the boundary from what one expects it 
to be there. In the present approach this problem can be avoided because the 
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boundary is entirely outside the physical space-time. From the causal properties 
of the evolution equations it is at least plausible that null infinity J which is 
a characteristic surface for the differential equations also numerically acts as a 
barrier for perturbations generated in the unphysical space-time. 

The section 4 is devoted to a discussion of the various gauge source functions. 
This is an important subject also in the conventional treatments of the Einstein 
equations because it is not clear what the implications are if one chooses e.g. the 
lapse function and the shift vector in that way and not in another. The existence 
of J imposes some questions concerning the resolution of features inside the 
physical space-time during the course of the computation. But those questions 
can be solved completely satisfactorily by making use of a certain choice of shift 
vector which is forced by the structure of conformal field equations. 

Finally, in section 5 we discuss the procedure of "radiation extraction" . By 
this term we mean the determination of certain asymptotic quantities which in 
a well defined sense characterise the gravitational radiation generated inside the 
physical space-time and escaping out to infinity. Here, this is a well defined 
procedure which involves finding the zero-set of the conformal factor, the in- 
terpolation of the field variables and a frame transformation to a well defined 
frame which is adapted to the geometry of J7. Another issue discussed in this 
section is the determination of the Bondi mass. Here, due to the different global 
topology, the situation is different from the physical one in that one cannot find 
a Bondi four-momentum but only a Bondi scalar. 

We tested the code by using exact solutions to first provide initial data for 
all the variables and then to compare the computed variables with the ana- 
lytic ones. We also computed the radiative quantities from the exact solutions 
for comparison with the numerically determined ones to check the radiation 
extraction method. 

As in [1] we use the conventions of [3] throughout. 

2 The numerical evolution scheme 

The evolution part of the conformal field equations presented in [1, section 
4] with the symmetry reduction described in [1, section 5] are to be solved 
numerically. To this end we set up a two-dimensional grid with coordinates 
u, V. It is assumed that the fields are periodic in u with period 2 so that we 
need to impose periodic boundary conditions at the surfaces u = ±1. The 
coordinate v takes values in the interval [—vo,vo]- The question of boundary 
conditions at the surfaces v = ±vo is rather delicate from the numerical point 
of view and will be discussed in detail in section 3. 

Having set up the grid we need to obtain solutions of the constraint equa- 
tions in order to initialise the fields. This should be done by solving the con- 
straint equations numerically given the appropriate free data and boundary 
data. However, since we are here concerned mainly with the evolution algo- 
rithm wc will simply take these initial values from appropriately adapted exact 
solutions which have recently been pointed out by Schmidt [2]. In particular, 
we take the rescaled A3 solution and one other solution from this class (see 
appendix A) in various gauges as our test case. 

As our numerical method for solving the evolution equations we choose finite 
difference schemes which are second order accurate in both time and space. In 
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particular, wc use the leapfrog and the Lax-WendrofF schemes both extended 
in a straightforward manner to two space dimensions. Of course, various other 
methods could have been employed. It soon becomes apparent that the leapfrog 
method is not a viable choice in our case. Since the conformal field equations 
form a quasi-linear symmetric hyperbolic system it follows that the character- 
istics which determine the evolution of the fields depend on the solution itself. 
Or, to put it differently, the wave parts of the fields propagate along the light 
cone which is determined by the metric which, in turn, is evolved by the field 
equations. Since the methods we employ are explicit we need to make sure 
that they remain stable by controlling the size of the time step At during the 
evolution. However, changing the time step dynamically cannot be done with 
the leapfrog method without loosing the second order accuracy in time. Hence, 
we will focus here exclusively on the Lax-WendrofF method. 

The equations in the general form given in [1, section 4] are manipulated 
using the NPspinor package [4] of Maple, extended to include the space spinor 
formalism. The equations are expanded into components using the decompo- 
sition into irreducibles for each spinor field. We should point out that the 
equations when written in components turn into, in general, complex equations 
for complex variables. Due to the reality properties of the spinor fields these 
equations come either in complex conjugate pairs or as real equations. This 
fact reduces the number and the complexity of the equations. The symmetry 
conditions given in [1, section 5] are used to simplify them. Maple is also used 
to test the equations thus obtained quite extensively in various ways: 

• They are checked against hand calculations in simple enough cases. 

• Inserting exact solutions into the equations should result in identities. 
These solutions are obtained also with the help of Maple using differ- 
ent routines by conformally transforming simple vacuum solutions of Ein- 
stein's equations with arbitrary conformal factors. 

• The most important test, however, is the fact that the evolution equations 
have to propagate the constraints. This property was verified for the full 
expanded evolution equations and constraints using Maple. 

The two-dimensional Lax-Wendroff scheme is a straightforward generalisa- 
tion of the one dimensional case. For a quasi-linear equation of the form 



it can be characterised as follows. Define the following operators acting on a 
grid function fij 



Diifkj = {fi+ij - fi-i^,,) , D,[f],j = (/,,^.+ . - . (2.3) 



ft = A{t, u, V, f)fu + B{t, u, V, f)fy + E{t, u, V, f) 



(2.1) 




(2.2) 
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Then the 2D-Lax-Wendroff scheme consists of the four steps 
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The further gcneraUsation 



of this scheme to three dimensions is straightforward. However, it becomes 
rather inefficient and it is here where probably operator splitting methods should 
be used. The complete discretisation of the equations was also carried out 
symbolically. 

The exact solutions described in appendix A have been used to provide the 
initial data and (in some cases) those boundary data which can be specified 
freely. It is clear from the form of the metric that these solutions have two 
space-like Killing vectors. The Killing vector dy is taken to be the one that 
is factored out by the symmetry reduction (see [1, section 5]). The metric 
functions are independent of x and y. If we choose z and x to correspond to the 
two remaining coordinates v and u respectively then the code is essentially a 
one-dimensional code. In order to test the two-dimensional performance of the 
code we therefore have to "warp" the coordinates. Thus, we put 



X = u, z = V — a{vQ — v^) sin(7ru). 



(2.4) 



We choose a in such a way that this transformation is bijective in the range 7i €E 
[— 1, l],v G [—vq, Vq\. In this coordinate system the orbits of the second Killing 
vector are distorted and not aligned with the grid, see Fig. 1. All computations 
which are presented here have been performed with these warped coordinates. 
We will now describe the criterion to determine the maximal time step At 
possible to evolve from an initial time level Sq at time t = t„ to the next 
time level Si at time t = tn + At. This is not specific to the Lax-Wendroff 
method but can be used with any explicit evolution scheme. The Courant- 
Priedrichs-Lewy condition [5] can be phrased as stating that "the numerical 
domain of dependence should enclose the analytical domain of dependence" . 
Now consider a point P in the future time level (cf. fig. 2). Its numerical domain 
of dependence consists of the points at time tn which are used to compute the 
field values at P. They lie within a rectangular area TZ bounded by coordinate 
lines u = u± and v = v±. The analytic domain of dependence is given by the 
intersection of the backward light cone of P with the time slice <So . The maximal 
allowed time-step At is therefore at most so big that the light cone just touches 
the boundary of TZ. To obtain a formula for the maximal At we analyse this 
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situation to first order or, what is the same thing, in Minkowski space. Then 
the time levels are planes and the light cone is a true cone, the null cone of the 
point P. Let O be the point in So which is "straight below" P in the sense that 
UP is a multiple of the normal vector of So- Wc take O as the origin. Lot Q 
be the point in So with the same spatial coordinates as P so that QP = Atdt. 
The equation for the plane <Si is given by {dt, x) = {dt, QP) = At so that 
Qpa ^ TvrAt. Then O^ is proportional to the shift vector O^ = -NAtT'di. 
The equation for the null cone of P is 

g^b {x" - OP'') (x^ - OP'') = (2.5) 

and its intersection E with <So is given by all points a;" which obey the equations 

taX" = 0, x"xa = 2N'^At^ = 0. (2.6) 

Now we take any plane H orthogonal to Sq, whose equation is LOaX'^ — s for 
some CO with Wa^" — and arbitrary s. Suppose that H touches E in a, point 
X. Then, at X the following equations hold: 

= of", X°-Xa = -2N'^At^, UJaX"- = S, Xa = Q!Wa- (2.7) 

The last equation expresses the fact that the tangent plane at X to in <So is 
parallel to the intersection of H with .Sq. From these equations, one can easily 
derive the relations 

2N'^At^ujaUj'' = -s^ (2.8) 

a;° = --N^At'^u". (2.9) 
s 

We are interested in coordinate planes within Sq. These arc obtained from uj — 
dx^ — NT^dt. In particular, we consider coordinate planes which are a distance 
±Ax^ from the point Q. These satisfy the equation Uax'^ = —NT^At ± Aa;*. 
Inserting this value for s in the equations above we obtain the equation 

(T* ± V-2a;aC^") NAt = ±Ax\ (2.10) 

which holds whenever the past light cone of P touches a coordinate plane which 
is ±A.t' from the point Q. Thus, according to our criterion, a valid time-step 
At should satisfy 

Ax* 

Ai<min— — , (2.11) 

where the minimum is taken over all points in So- There are some points to be 
mentioned: 



• In our present case the square root is simply y/—2uJaOJ°' = 2^fC\^fi^^ , so 
that determining the maximal time-step is rather simple. It involves going 
through the grid and finding the maximum of some algebraic function of 
the fields (no inversion of the spatial metric). 

• We find that the criterion is at least necessary for stability, i.e., if we do 

not enforce the time-step to be at most the above value then the scheme 
becomes unstable. So far it has also been sufficient for stability. 
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• This criterion might be eonservative. In fact, one could imagine that one 
should be able to increase the time step until the first of the adjacent grid 
points comes to lie on the null cone, the others still being not inside. We 
have not investigated this further. 

• This is a first order criterion and it might be too crude for the Lax- 
Wendroff method. One could think of enforcing this condition at each 
half step. Again this has not been investigated. 

A complete time step is performed by going through the following steps given 
the solution at the time level f„ 

• Find the maximal possible time step At by inspection of the current time 
level. 

• Set the gauge functions, also possibly according to the properties of the 
solution at the current time level. 

• Update the solution at the interior points using the gauge functions and 
the time-step by performing the above four steps for each function. After 
the first half step specify the gauge source functions again. 

• Update the solution at the boundary points. 

3 Boundary treatment 

Analytically, the hyperboloidal initial value problem does not need any bound- 
ary conditions. The initial data arc given on a three dimensional manifold S 
with boundary dS on which the conformal factor Q is supposed to vanish. Then 
a solution exists on the four dimensional manifold M = 5 x [0, t] for some r > 
which is such that the boundary dM = dS x [0, r] is a null hypcrsurfacc and 
hence characteristic. That means that even if one would extend the initial data 
across the boundary f2 = in some way this extension could not influence the 
interior, i.e., the physical space-time depends only on the data given inside the 
boundary. 

The situation is different in the numerical case. The characteristic speeds 
arc different for different modes which arc propagated by the numerical scheme. 
In particular, non-physical modes tend to propagate at speeds much higher than 
physical propagation speeds and thus contaminate the solution all over the en- 
tire computational domain. A notorious place where non-physical modes are 
generated is at the boundaries of the domain. Due to the lack of enough grid 
points there, in general, the numerical evolution scheme has to be changed. It is 
absolutely vital to impose boundary conditions so that the non-physical modes 
are kept small. The GKS theory [6, 7] which has been developed for analysing 
such situations is inherently difficult to apply. A different (equivalent) formu- 
lation based on the notion of group velocity for finite difference schemes has 
been given by Trefethen [8]. It has been found that certain intuitive numerical 
boundary conditions do not perform as expected. Conditions which work for 
one numerical scheme do not necessarily work for others. For linear equations 
in one space dimension the mathematical analysis can completely be carried 
through. It turns out that the essential criterion is a non-degeneracy condition 
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for a linear system of equations obtained from the combination of the evolution 
scheme and the boundary condition. This system is required to have no solu- 
tions in order to exclude the parasitic modes. Although Trefethen's method is 
very physical and intuitive it does not provide enough information in the case of 
higher dimensional and/or non- linear equations. It does, however, give valuable 
hints as to which conditions might have a chance to be useful in those more 
general situations treated with the Lax-Wendroff scheme. 

The situation is somewhat ironic in the present case. One is not at all 
interested in what happens at the boundary because this is (usually) outside 
the physical space-time. However, it is the boundary which needs the most 
careful treatment. One would wish to find gauge conditions which make the 
non-physical portion of the computational domain small, ideally putting the 
boundary at f7 = 0. We will examine the feasibility of this idea later on. 

In another aspect, the present situation is also quite disadvantageous. Usu- 
ally, it is of great importance that the analytical problem has a well posed 
initial-boundary value problem. The rigorous analysis provides the information 
about which data can be specified freely at the boundary and which data is de- 
termined from information propagating towards the boundary from the inside. 
Knowledge of this kind is necessary in order not to over-specify the solution 
at the boundary, because this would inevitably lead to instabilities (except for 
extremely simple cases). In our case, it is not known at present whether the 
system admits a well posed initial-boundary value problem (see, however [9]). 
To overcome this lack of information we analyse the system to first order at the 
boundary in the following sense. 

The boundary ,B is a time-like three-dimensional hypcrsurface in the space- 
time. Let Ua be the space-like conormal oi B. A system of partial differential 
equations which has the form 

dtUa + A'^i^diu'^ = ba (3.1) 

can be rewritten in the form 

dtUa + C^aVu^ + B^f^dAU^ = ba, (3.2) 

wh(^r(^ V is the derivative along any vector field i//' which satisfies Uau"' = 1 on B 
(usually taken to be the normal vector field extended off B in an arbitrary way) 
and where the Oa are derivatives intrinsic to B at points of B. On the bound- 
ary the matrix Caf} regulates to first order the propagation of the fields across 
the boundary. By analysing its structure one can gain valuable insights into 
the behaviour of the solution on B. In particular, finding the eigenvalues and 
eigenvectors of Caf3 (which in our case is hcrmitian) enables us to select com- 
binations of the fields which (to first order) propagate purely inward or purely 
outward or which stay on the boundary. These have to be treated differently. 
While the ingoing pieces can be prescribed freely, the outgoing ones have to 
be obtained from the interior. This is done here by extrapolation. That this 
might be possible is indicated be the Trefethen analysis which shows that ex- 
trapolation remains stable when used in conjimction with the one-dimensional 
Lax-Wendroff method. We want to mention that this analysis applies not only 
at the boundary but also at the interfaces between grid cells. This is impor- 
tant for possible future application of high resolution methods which require the 
solution of Riemann problems at each grid cell, see e.g., [10]. 
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To be somewhat more precise we need to analyse the three subsystems of 
the full system which do not consist entirely of advection equations along the dt 
vector. These are the systems for the variables {Kabcd, Kab, K), the variables 
{4iABCD,(t>AB,4>) and for the Weyl curvature {Eabcd, Babcd), respectively. 
Note, that this analysis is valid in the three-dimensional case. Only in the 
code we have specialized this to the two-dimensional case. Let us describe the 
procedure for the 0-system. 

First, we note that Ha = uab can be viewed as a complex metric on spin 
space which reduces the structure group from SL{2, C) down to U{1). Thus, it 
is possible to decompose any symmetric spinors $ab, ^abcd into irreducible 
pieces with respect to the smaller structure group: 

^AB = ^%-^,nAB^^'\ (3.3) 
1 3 

^ABCD = '^%CD - ■^HAB'^CD) + ^"(AB^lcD) (3-4) 

where UABn^^ = —2r? and 

$(1) = $Asn^^, (3.5) 
$^'^=$AscDn^^n^^, (3.6) 

^% = ^ABCDrf'' -r^ncD^^^y (3.7) 

This decomposition which is very natural algebraically corresponds geometri- 
cally to a decomposition of the fields into parts which are vertical and tangen- 
tial to B. The principal part of the ^system which does not contain tangential 
derivatives is 

d(f>ABCD - n(^AB'D<t>CD), (3.8) 
d(l>AB - yiAB'Dcj) + n^'^V(t>ABCD, (3.9) 

d(j) + 2nAB'D(t>'^^. (3.10) 

Here we have neglected terms containing derivatives of uab because those do 
not alter the symbol of the subsystem. Inserting the decompositions of the 
variables we get the following system 

8(1) + 2V<j>'-^\ (3.11) 
a0(i) + l„2p^ + 2:)<^(2)^ (3.12) 

O 

a0(2) + 4 2^^(l)^ (3.13) 
o 

8cj)f^+V<p% (3.14) 

8cj)^^l+n^Vcl>f^, (3.15) 

dcl^%cD- (3-16) 

Obviously, this system can be decomposed into three smaller subsystems and 
it can be shown that the coefficient matrix of the V operator is hermitian with 
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respect to a suitable inner product (it has to be because it comes from a sym- 
metric hyperbolic system). Now it is easy to find "characteristic combinations" 
of the variables so that the symbol becomes diagonal, i.e., it has the form 

dC + XVC (3.17) 

for each characteristic quantity. These combinations are unique only up to a 
scaling factor. We choose the following quantities with their respective charac- 
teristic speeds 

Co = 4nV - Q(t>^'^\ A = (3.18) 

C± = 4nV ± 3\/4^(?!)^^^ + 3(A^^\ A = ±\/4^ (3.19) 



C±s = ± V^^(fil A = ±Vn2 (3.20) 

Cabcd = cP%CD^ A = 0. (3.21) 

In an analogous way we find the characteristic quantities for the JiT-system 

Co = 3^(2) + 4^2^^ ^ ^ Q 22) 

C± = 6K^2^ ± sVI^k'-^^ + 4n^K, A = ±V4^ (3.23) 



CIb = ^K^AB ± ^^n^Z, A = ±V2n2 (3.24) 



^AB ~ ^^^AB 

Cabcd = K-fj^^D^ A = 0. (3.25) 

The Weyl system has a completely different structure from the previous two 
subsystems. Nevertheless, the analysis yields characteristic quantities written 
in terms of the complex Weyl spinor tpABCD and the spinor field i^abcd = 

n(A^1pBCD)E- 

Co=V'°\ A = (3.26) 

C±5 = ±x/^Vi'i-<^«, X = ±V^ (3.27) 

^± _ _L,/:;2./,(o) o, JO) 



CIbcd = ±^n^rih - 2V^% , A = ± V4n2. (3.28) 



In our case, the boundary is given as a surface v = const so that we can put 
HAB = C\g. Let us now focus on (3.17). Inserting the explicit expressions for 
the derivatives this is 

dtC + N {-T'^dyC + XdyC) (3.29) 

Therefore, it is the sign of (A — T^) which regulates in which direction the 
quantity C propagates across the boundary. 

To update the values at the boundary points we proceed as follows. First we 
determine the characteristic quantities on the boundary. This is done by looking 
at the sign of the corresponding eigenvalues which decides whether to simply 
set the value arbitrarily in case the quantity propagates inwards or else whether 
to find the value by extrapolation from the interior. Prom the characteristic 
quantities we obtain the field values by reversing the transformations above. 

In the situations considered this procedure yields a stable algorithm. This is 
consistent with the Trefethen theory which shows that in the one-dimensional 
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case extrapolation together with the Lax-WendrofF time evolution scheme re- 
mains stable. By its very nature our procedure is a first order approximation 
to the real situation so that we cannot expect to obtain a code which is sec- 
ond order accurate in the neighbourhood of the boundary. However, since the 
surface f2 = is a characteristic surface we may hope, that the error does not 
too severely influence the physical space-time as long as the boundary of the 
computational domain is outside the physical space-time. 

4 Gauge choices 

In this section wc want to present some results about the various possible gauge 
choices. Our emphasis will be on the temporal gauge choices. There is one 
class of choices for the shift vector which is natural in the present context of the 
conformal field equations. The gauge for the frame rotations and the third class 
of gauges, namely the choice of the scalar curvature A, is unknown territory as 
of yet and we put the corresponding gauge source functions equal to zero in the 
code. As was pointed out in I, in the; c;ase of the frame rotations this means 
that the frame is Fermi- Walker transported along the normal vector of the time 
foliation. 

4.1 Choices of lapse 

Let us start with the temporal gauges. Fixing the lapse function is a difficult 
task. This function has to be chosen in such a way that the time coordinate 
does not degenerate in the course of the evolution. Here is an attempt to collect 
some criteria which should be satisfied by the lapse: 

• the lapse function should not "collapse" in the sense that it approaches 

zero in a finite coordinate time, 

• the surfaces of constant time should remain smooth, 

• the lines of constant spatial coordinates should not intersect, 

• the lapse function should remain positive, 

• it should not develop too steep gradients, 

• depending on the problem the foliation should or should not avoid singu- 
larities. 

In our treatment of the hypcrboloidal initial value problem, the lapse function 
cannot be chosen directly. Instead it is governed by an evolution equation which 
contains the "harmonicity" F = 2Dt as an arbitrary function of the coordinates. 
It is not easy to find a function F which allows the lapse to satisfy some or all 
of the above criteria. The reason is that one has no idea what F should look 
like in coordinates which are constructed as the code moves along. To make the 
lapse satisfy the criteria above one needs a certain amount of "feedback" i.e., 
information about the current status of the evolution seems to be unavoidable. 
This means, that one should specify as a function of the field variables. But 
since in the system also the derivatives of F appear this leads to the problem that 
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the characteristics of the system change because the symbol has been altered. 
We will discuss this later in this section. 

The various choices for N that have been considered so far are 

• the "natural gauge" for the exact solutions obtained by setting F equal to 
the expression computed from the explicit form of the metric (cf. appendix 
A) and variations thereof, 

• the "GauB gauge" which is the condition that N should be constant 

throughout the timcslice, 

• the harmonic gauge with F = and 

• the special class of gauges for which F is a function of N and K only, 
F = f{N, K), which in fact includes all of the above gauges. 

The popular "maximal gauge" where one requires K = probably cannot be 
achieved by specifying F. 

We can judge the effects of these gauges by monitoring the function t which 
satisfies the eiconal equation VaTV°r = 1 and the condition r = on the initial 
surface. The value of this function at a point P gives the distance in proper 
time between P and the intersection point of the geodesic through P tangent to 
the unit normal of the foliation and the initial surface. This function is evolved 
simultaneously with the other field variables. For the A3 solution the proper 
time distance from a point on the initial surface t = to with z = and the 
singularity att = z = Oisr = 2^— to for to < and we can compare how 
far the evolution proceeds with the various gauge choices. In all the examples 
presented in this section a 100 x 100 grid was used and the coordinate system 
is one discussed in section 2 with tig = 5. The initial data wore taken from 
the exact A3 solution at an initial time to = —5. The boundary values were 
specified depending on the cases considered. When the "natural gauge" was 
used the ingoing boundary data were taken from the A3 solution. In all other 
cases these values were used initially to satisfy the corner conditions and then 
specified to decrease exponentially, so that after approximately 20 timesteps the 
ingoing values are zero. 

4.1.1 F as a function of the space-time coordinates 

For the "natural gauge" we find that we can in principle approach the singularity 
arbitrarily closely by increasing the resolution appropriately. Obviously, there 
are hard limits imposed on the calculation by the finite precision arithmetic so 
that eventually the numbers will overflow. 

The difference between the natural and the harmonic gauge is shown in Fig. 3 
and Fig. 4. In both figures we show the proper time r and the lapse at a late 
instant of the time evolution. For the natural gauge this was dictated by the 
code which stopped because the time step could not be chosen without violating 
the CFL condition. This is due to the closeness of the singularity. The light 
cones arc infinitely stretched along the symmetry directions as the singularity 
is approached. A better resolution would have extended the evolution time 
somewhat more but eventually the same thing would happen. In this run the 
coordinate time elapsed was roughly 4.405. It is clear that with this temporal 
gauge the singularity can be reached (at least in principle). 
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With harmonic gauge the code was stopped after an elapsed coordinate time 
20.07. In Fig. 4 it can be seen that the evolution close to the singularity has 
slowed down considerably. The proper time in the center is much smaller than 
for regions further outside. The lapse function is very small in the center and it is 
decreasing. It is not known whether the evolution will reach the singularity even 
in principle. The decrease in the lapse could be so rapid that the integrated 
proper time along the central geodesic would reach a limit below the value 
2^y—to. It is quite likely, that in this gauge the code will ultimately not be able 
to resolve the steep gradients which occur at late times between the interior 
region which cannot progress beyond the singularity and the exterior region 
which can. 

Related to this phenomenon is the fact that the interior region shrinks. On 
the initial surface is located at the boundary of the grid and it gradually 
moves inward during the evolution. Ultimately, there will be only very few grid 
points left in the interior region. This is a phenomenon which has nothing to do 
with the temporal gauge but with the choice of the shift vector. We will discuss 
later in this section a shift gauge which allows the freezing of J on the grid. 

To get some feeling for the influence of F on the time slicings we study the 
slicings obtained from substituting p ■ F for F with some parameter p. For p = 
we have the harmonic gauge which avoids the singularity. For p = 1 we have 
the natural gauge which allows us to reach the singularity in finite time. What 
happens when we increase p beyond unity? 

We evolve for a fixed coordinate time interval t G [—5,-4] with different 
values of p G [1.0, 1.8315]. We find that for p > 1.8315 the code crashes. It 
is easily seen that this crash cannot be due to the curvature singularity in 
the A3 solution but is a coordinate singularity. In Fig. 5 we plot the proper 
time distance between the initial and final time slice along the central geodesic 
(z = 0) versus the parameter p. We see that t(p) has an infinite derivative at 
p = 1.8315 with a finite value of t far from its value 2-\/5 at the singularity. 
The lapse function N for different values of p diverges rapidly (cp. Fig. 6). 
The curvature invariant I = fl^ i^EoEi + Gi?!) which diverges at the singularity 
stays perfectly regular. Fig. 7 shows the exact invariant plotted against proper 
time along the central geodesic. The dots are the values of / and r obtained 
from the runs with different parameter values. We see that the behaviour of 
these functions is not altered by the occurrence of the coordinate singularity. 

Finally, we show in Fig. 8 the profiles of the conformal factor for various 
values of p. As the parameter approaches its final value the conformal factor 
develops a minimum at the center. Although this behaviour seems strange 
at first sight it can easily be explained. The conformal factor decreases as a 
function of r for fixed spatial coordinates. Due to the rapid divergence of the 
lapse the proper time at z = is much larger than in the regions outside. So 
that we see values of O in the center which are reduced over proportion from 
the values outside the center. This accounts for that central dip. 

4.1.2 The Gaufi gauge 

The "Gaul3 gauge" which forces N to be constant is a condition which is im- 
posed on the lapse function directly. In principle, it is possible to express the 
exact solutions in Gaui3 coordinates by performing the coordinate transforma- 
tion explicitly. Then one can compute the "harmonicity" for these coordinates 
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and do the evolution. However, we proceed somewhat differently to impose the 
Gaufi gauge. The lapse function satisfies the evolution equation [1, (4.29)] 

dN = -N'^K - N^F (4.1) 

where d = Nd. Now suppose that A'' would satisfy an evolution equation 

ON = a{N - N), (4.2) 

where N is an arbitrary (positive) constant. This equation has the solution 
N{s) = N + TVpe""", s being the parameter with ds = 1. Thus, for a > the 
lapse N approaches the constant value N in the course of the evolution. We 
can make A'' satisfy the above evolution equation (4.2) by choosing 

F=^fc^±^. (4.3) 

We find what one would expect, namely that in this gauge the timeslices develop 
caustics (or, what has become known as "coordinate shocks"). This makes the 
Gaufi gauge inappropriate for long time evolutions. 



4.1.3 F as a function of N and K 

Let us now focus on the class of gauges defined by _F = f{N,K). Among 
this class there is a subclass for which the lapse depends only on the three- 
dimensional volume (-element) V ^ N = g{V). For these gauges, we have with 
appropriate assumptions on the function g 

dN = g'{V) dV = -Kg'iV) V = -Kg'{g-\N)) g-\N), (4.4) 

so that, in fact, is a function of and K only. The natural gauge falls into 
this subclass with g{x) = x~3 and, consequently, F = — Similarly, the 
harmonic gauge with F = is in this class with g{x) oc x. 

If we specify F for the natural gauge not as a function of the space-time 
variables but instead as depending on the field variables, then an interesting 
phenomenon occurs. Although nothing else in the code has been changed it 
seems to notice this difference because the boundary becomes unstable very 
quickly. However, inside the computational domain we get the same solution 
without any significant difference between the two ways of specifying the gauge. 

This phenomenon can be traced to the fact mentioned above, namely that 
the gauge specification might change the characteristics of the system. We can 
see this explicitly as follows. In equation [1, (4.31)] the derivative of F appears. 
With F = f{N, K) the principal part of that equation is 

OKab + 2d^'^KABCD + 2NlKdABK + 2NlNdABN. (4.5) 

The term involving OabN can be removed by using the constraint equation [1, 

(4.44)] so that the symbol for the X-subsystem is {p,Pab) Sp{L, K), where Sp 
denotes the sesquilinear form obtained from the principal part of the JC-system 
by replacing the derivative operators (OjOab) with {p,pab) and multiplying 
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appropriately with the complex conjugate of some spinor fields {Lab, Labcd)- 
Thus, 

Sp{L,K) = I L^^'Kab + L^'^P^^'Kabcd + {NIk)L^''pab K 

+ pL^^^'^Kabcd - L^^^'^pabKcd. (4.6) 

Various important properties of the -ff-system can be determined from the form 
Sp. In particular, the system is symmetric if Sp is hermitian, i.e., if Sp{L, K) = 
Sp{K, L) for all {p,pab)- It will be symmetric hyperbolic iff there exists {p,Pab) 
such that Sp is positive definite. We see from the above that the Jf-system 
can be made symmetric if we do not consider equation [1, (4.32)] but add to 
that equation an appropriate multiple of its trace. This changes Sp into (with 
L = L^^ab) 

Sp{L, K) = ^ L^^'Kab + 2^ Vi^ABCiJ + {Nf,K)L^''pAB K 

+ pL^^^''Kabcd - L^^^'^pabKcd 

+ {NIk)p LK - {NIk)L p^^Kcd. (4.7) 

It is easy to see that this form is hermitian and that it will also be positive 
definite provided that the inequality 

a = NfK + l/S>0 (4.8) 

is satisfied. This is of course a restriction on the possible gauges. 

Furthermore, the characteristics of the above system can be obtained by 
inspection of its characteristic polynomial defined by 

P{p,PAB)=det{Sp), (4.9) 

for which we obtain the expression 

P{p,Pab) = 2Aap^ (p^ + p^^p^s)^ (p'' + 2 (^a + ^ P^^'pab^ ■ (4.10) 

This polynomial is homogeneous of degree nine in its variables and, regarded as a 
polynomial in p only, it will have nine real zeroes provided that a + 1 is positive, 
which is always the case if the inequality (4.8) is satisfied. In this case, there will 
exist three different characteristics, namely the lines along the time evolution 
vector, the cone given by the first factor in parenthesis in (4.10) which is double 
layered and a simply layered cone given by the last factor in (4.10). The latter 
cone is gauge dependent while the former is not. The degenerate characteristic 
is time-like while the gauge dependent characteristic has no gauge independent 
causal character. The cases when F vanishes, the harmonic gauge, or when F 
is specified as a space-time function correspond to a = 1/3 in which case the 
gauge dependent characteristic coincides with the light cone. For the other cases 
< a < 1/3 and 1/3 < a the characteristic is time-like, respectively space-like. 
However, there are gauges within the specified class for which this characteristic 
does not even exist. Thus, the system acquires a mixed type, having hyperbolic 
and elliptic parts. In particular, for the natural gauge specified in terms of field 
variables we have Nf^x + 1/3 = —1 which violates the inequality (4.8). 
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A natural question to ask is the following: to what extent are these features 
noticeable in the code? Judging from experiments what seems to be the case is 
that the code will probably not detect differences in the various cases as long as 
it does not make use of the hyperbolic character of the system. In particular, it 
will probably not detect when the system changes its character from a hyperbolic 
type to a mixed type due to a gauge change. However, in those instances where 
the hyperbolic character is in fact used in the code difficulties will arise. In the 
present code we find that the boundary will become unstable very quickly when 
we choose a gauge which makes the system partly elliptic (we use this term 
only to indicate that the resulting system is no longer hyperbolic). This is of 
course due to our treatment of the boundary which implicitly assumes that F 
is specified as a coordinate function. Another instance which can detect gauge 
changes is due to the time-step control. Here, we implicitly assume that the 
largest propagation speed is the speed of light. For gauges with a > 1/3 this 
is not the case, the largest propagation speed is bigger than the speed of light. 
But the largest speed is the one which limits the time-step in order to enforce 
the Courant-Friedrichs-Lewy condition for stability of the code. And, in fact, 
choosing a big enough results in numerical instabilities inside the computational 
domain. 

These tests have been performed using the initial data of the A3 solution and 
then specifying various gauges by choosing F. We find the surprising feature 
indicated already above that the code detects whether F is specified as a space- 
time function 

F=^^= (4.11) 

or as a function of lapse and mean curvature F = — | ^ • While it runs without 
problems in the former case all the way up to a maximum of the proper time r 
close to its theoretical limit in the latter case the boundary becomes unstable 
very quickly. As surprising as this might seem it is still in accordance with the 
general picture. What might be even more surprising is the fact that in the 
interior there is apparently no sign of any difference between the two cases. 

Another gauge which has some geometric significance is given by choosing 
N oc \/V. This condition can be obtained from the requirement that the height 
of the backward light cone of a point in the next time level should be proportional 
to the "volume radius" R = ^/V of its intersection with the current time level. 
This condition is satisfied for the standard t-foliation in Minkowski space. Thus, 
we have 

(4 12) 

iV " 3 T/ " 3 ^ ' 

and F = I ^ . The speed of the gauge modes is in this case bigger than the 
speed of light but the system remains hyperbolic. In practice, this gauge is not 
very much different from the harmonic gauge. 

What wc learn from these various discussions and experiments is that the 
natural gauge is the most efficient one for approaching the singularity. However, 
in situations where there is no exact solution this gauge is not available. Now 
one has various possibilities: one could prescribe a gauge condition once and 
for all like the ones considered in section 4.1.3 or even like the maximal gauge 
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where one needs to solve an elliptic equation on each timeslice. The former have 
the disadvantage that they introduce superluminal propagation speeds into the 
problem so that the stability of the (explicit) methods forces rather small time 
steps while the latter arc rather time consuming. The other approach would 
be to always specify as a function of the coordinates. This means that one 
needs to experiment in order to find a good candidate expression for F which 
allows to reach singularities effectively. This method is very flexible but it is 
also rather obscure because there are no guiding principles about the shape of 
the harmonicity function F. 

4.2 Choices of shift vector 

The choice of a shift vector is even more obscure. There are two issues involved 
in the choice of the shift vector: the problem of what to do at the points of the 
physical space-time and how one is to treat the points on J. 

Let us first discuss the interior issues. To describe the problems involved we 
focus on the lines of constant spatial coordinates parametrised by the coordinate 
time. These are the integral curves of the vector field t = ^ , the "t- lines" , which 
form a family of time-like lines. It is the geometry of that congruence which 
can be influenced by choosing the shift vector. To discuss this in more detail 
we decompose the time-like coordinate vector into lapse and shift 

t" = N {t"- + T") (4.13) 

and we choose a connecting vector ^"'j i.e., a vector field which commutes with 
r". Such a connecting vector, which is also called a Jacobi field, can be viewed 
as describing an infinitesimally separated line in the family with f connecting 
points with the same value of the time parameter. Thus, ^" is tangent to the 
t = const surfaces, satisfying ^"'ta = 0. Prom the commutator of the two vector 
fields we obtain 

T^Wa^'' = ^rVaAfr'' + N (rV„i^ + rVaT*-) . (4.14) 

The contraction of this equation with the time-like normal of the surfaces yields 
the constraint equation [1, (4.44)] which couples the gradient of N to the time 
evolution of the acceleration vector. The other part of the cqiiation which is 
intrinsic to the hypersurfaces can be obtained by projecting (4.14) along r" onto 
the hypersurfaces. This is achieved by contraction with the projection operator 

P% = - -^^bT^ (4.15) 

This yields the relation 

C=N {^''Kb" + i^dbT'') , (4.16) 

where the dot simply means r^Va followed by projection. As in the case of 

geodesic congruences this family of f-lines can be described infinitesimally by 
its twist, shear and divergence according to the irreducible decomposition of 
the right hand side of (4.16). Prom this we can conclude that a constant shift 
vector generally causes the family of t-lines to shear and diverge, depending 
on the properties of the extrinsic curvature. This is well known in the case 
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of Gaufi coordinates which develop conjugate points unless the hypersurface is 
very special. 

The goal of choosing a shift vector should be to prevent the t-lines from 
coming too close together. The twist of the congruence, entirely due to the shift 
vector, does not change the relative distances of the f-lines. Therefore, we need 
to eliminate as many components as is possible from the shear and divergence 
combined in 



Since there are only three components in the shift vector, only three compo- 
nents of (Tab can be compensated. Depending on which components are to be 
eliminated there result different, and in general elliptic, equations to be satisfied 
by T". One possibility is to eliminate the divergence of the congruence which 
leads essentially to a Poisson equation. Another possibility to determine a shift 
vector is not to eliminate components of aab but to minimise the functional 
/ aabcr°'''dV. This leads to the well-known "minimal distortion" shift condition, 
which is a second order elliptic equation for the shift vector. The problems re- 
lated to the interior of J, i.e., to the physical space- time are essentially the same 
as in the numerical treatment of the traditional Cauchy problem, and there is 
no insight to be gained from the hyperboloidal initial value problem. 

However, this is difi^erent when one looks at the issues concerned with the 
boundary of the physical space-time. One objection against the use of conformal 
methods in the numerical treatment of the Einstein equations has been the 
following: as the evolution proceeds the part of <S which corresponds to the 
physical space-time shrinks so that there are less and less grid points left in 
the interior of J7 (see the figures 3 and 4). This implies that the resolution 
of features in the physical space-time is getting smaller. However, as it turns 
out, this is a misconception which might be caused by the familiar conformal 
diagrams of asymptotically fiat space-times. There it is assumed that light rays 
are aligned on 45° lines. This need not be the case. In fact, by choosing the 
shift vector on J7 appropriately we gain complete control over the movement of 
J through the grid. Therefore, we get to choose between (at least) two options. 
On the one hand, we can compute a Penrose diagram of the space-time which 
is useful for discussing its global properties. E.g., it helps in deciding whether 
there exists a regular i+ or whether there appear singularities before i+ can be 
reached. Another option is to have J not move at all through the grid. This 
enables one to keep the resolution in the interior constant so that the physical 
space-time docs not suffer any loss of resolution (luring the evolution. This 
property is desirable when studying the behaviour of sources in the physical 
space-time. Although in this case, the picture which emerges looks like the one 
obtained by spatially compactifying space-time one should keep in mind that 
the conformal structures are entirely different in the two cases. After all, in the 
picture proposed here, J is still a regular characteristic surface. 

How can we achieve that does not move through the grid? The equation 
for the conformal factor is 



(4.17) 



dtSl = N {T'difl + S) . 



(4.18) 
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Note, that T%n = T^^Bab^ = T'^'^'Sab- Thus, if we choose 

^AB 



Tab = (4.19) 



then we obtain the equation 



d,n=^(^EAB^^'' + l^') (4.20) 

4nN 



(5-f2A). (4.21) 

Therefore, 9(fi is proportional to fl so that fl remains zero along the t-lines at 
those places where it was zero in the beginning of the evohition. i.e.. on J'. This 
imphes that J does not move through the grid. Although it looks as if the shift 
vector is now uniquely fixed, this is not the case. Note, that the choice 

Tab = 2^ + ^Tab (4.22) 

exhibits the same behaviour. Here Tab is completely arbitrary apart from the 
fact that it should be bounded on J'. Its only effect is on the coordinates in 
the interior. If we choose Tab so that ^Tab has finite values on J then we 
can achieve that J moves through the grid in a rather arbitrary but controlled 
fashion. 

It should also be pointed out that the form of the shift vector given in (4.22) 
is unique, imposed by the geometry. It does not sufii'er from the shortcomings 
of other gauge choices. In fact, although it is specified by prescribing Tab as 
a function of the dependent variables, this does not change the characteristics 
of the system even though there are terms involving the derivatives of the shift 
vector. 

In Fig. 9 we show the proper time and the lapse function for a run with 
harmonic gauge and scri freezing. The initial location of JT" was on the boundary 
Vo = ±5. The length of coordinate time spent was ti — to = 25 with roughly 1000 
time steps. has moved during this evolution at most over 3 grid points. This 
is due to numerical inaccuracies. We see from the figures that the evolution is 
much more homogeneous over the interior with differences in proper time within 
the interval [2.96,3.04]. But we also see that the lapse has decreased rapidly, 
from a maximum value of 0.316 at the beginning to a maximum value of 0.0422 
in the end. 



5 Mass loss and radiation 

The main motivation to consider the conformal field equations in the first place 
is the claim that having J' at finite places allows a well defined numerical de- 
scription of the asymptotic properties like the radiative information (such as 
shear and news on J') and also the global properties like the Bondi energy- 
momentum and angular momentum. From the nature of the hyperboloidal 
initial value problem it is clear that we cannot get our hands on the ADM 
quantities which are located at space-like infinity z° which is not in the domain 
of dependence of any hyperboloidal initial surface. 
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In the numerical treatment there exists a natural foliation of J7 into two- 
dimensional cross sections or "cuts" which is obtained from the intersection of 
J with the constant time hypersurfaces. The news, shear and the "null datum" 
are local quantities in the sense that their value at a point on is constructed 
from the values of the field variables at that point. Therefore, these quantities 
are not sensitive to the topology of the "cuts" . In contrast, the Bondi quantities 
are global concepts and there is currently no way to determine their value from 
only local information. As a consequence they are very sensitive to the topology 
of the cuts and, in fact, they are so sensitive that in our case study with "cuts" 
which have a toroidal topology there does not exist an energy-momentum four- 
vector but only one scalar quantity which we still call the Bondi mass. The 
reason behind this unexpected phenomenon will be discussed below. 

The main problem one is faced with when trying to obtain expressions for 
the asymptotic quantities is the fact that J does not look the way it does in the 
analytical treatments. In particular, there one usually assumes that a conformal 
gauge has been chosen so that J is divergence free, i.e., that the area of a cut 
does not change when the cut is moved along the null generators of J. Since 
J itself is shear free this implies that the shape of a cut does not change either 
along and this fact can be exploited to choose the metric on that family of 
cuts to be one with constant curvature. Usually this is the unit sphere metric. 
In our case, where the cuts have toroidal topology one would choose a flat metric 
on the cuts. 

In a numerical treatment where the conformal factor is one of the evolving 
variables one has almost no control of its behaviour (at least at present). Thus, 
we do not have the freedom to specify that should have these nice properties 
and have to live with the way it emerges from the numerical computation. The 
only way to possibly influence the behaviour of the conformal factor is by way 
of tuning the gauge source function A for the conformal gauge. However, this is 
a rather indirect way and at the moment it is completely unclear whether (and 
how) one should specify A so that J does have the desired properties. 

Another point is that the radiative quantities are referred to a specific tetrad 
(or spin frame) on J' which is adapted to the geometry there. Again, in the 
numerical treatment the tetrad is fixed by other means which implies that we 
need to transform from the given tetrad to the geometric tetrad in order to 
obtain the correct values of the asymptotic quantities. Again, it is not yet known 
how to impose gauge conditions so that the computed tetrad always coincides 
with the geometric tetrad on J^. The transformation from the numerical frame 
to the adapted frame is straightforward. Recall that the condition imposed on 
the adapted spin frame {0^,I^) is [3]: 

Vaa'O = -AIaIa' + 0{Q). (5.1) 

This condition fixes the direction of the null vector I"^!^ but says nothing 
about the space-like vector 0^1^ and its complex conjugate. Given a cut of J 
these are required to be tangent vectors to the cut. Then the transformed spin 
frame is fixed up to the scalings O"^ i— > c- O"^ , 7"^ i— > • 7"^ . The transformation 
to the new spin frame is 

= ao^ + bL^, (5.2) 

I^ = ^^A-bo^ + ai% (5.3) 
aa + 00 
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with a = —c^^p- and b — c (^^^ — l) and c an arbitrary complex function on 
the cut. With this choice of a and b we have achieved that on J the null vector 
I^I"^ is aligned with the null generators of J. The factor A in (5.1) is found 
to be 

A = 4cc (^^21 (5.4) 

and we furthermore fix c = ^ for the remainder of this section. 

The asymptotic quantities with respect to the adapted frame can now be 
expressed on J in terms of the field variables. These expressions are rather 
lengthy in the general case but quite manageable in the symmetry reduced 
case that we are looking at here. Following is a list of the variables which are of 
interest to us and the expressions to compute them in terms of the field variables 
in the reduced case: 

cr' = 0, k' = 0, (5.5) 

P =-^, r (5.6) 

(7 = -I (K40st2 + K44 + 6K42SI2) (5.7) 
o 

V'2 = J {E0SI2 - 2E2 + E4SI0) , (5.8) 

V'4 = E4sio + 4 Baslo + 6 Eisl^ + 4 B1S20 + Eq, (5.9) 

=2 ('?^44s|o + 6(?i42sio + (f>4o) (5.10) 

where we have introduced S22 = and S20 = • The function Af is the so 
called "news" function. 

Having these expressions at hand it is in principle straightforward to obtain 
the asymptotic quantities from the numerical data. The only obstacle is that 
the level sets of fl do not necessarily agree with grid lines so that one has to 
trace out the zero set of ft within the grid and then interpolate for the values 
of the field variables there. This task is greatly simplified by using the J fixing 
shift gauge discussed in section 4 when it is possible to align J on a grid-line 
initially. 

In Fig. 10 we present a surface representation of the null datum tp4 for the A3 
solution. This is a non-radiating solution so tp4 should vanish. Indeed, we find 
that only when the singularity is approached the function differs significantly 
from zero. This is due to the closeness of the singularity. We should also 
point out that this figure has been produced in the warped coordinate system 
without the use of the J fixing shift gauge. It is only in the late stages and 
in the central region where the warping is maximal when the tracing out of J 
produces too large errors. In a similar way the Wl solution was treated. Now, 
V'4 cannot be expected to vanish because this is a radiating solution. Since there 
is an additional symmetry present in the solution which is aligned along J the 
function should be constant on J'. We found that the tracing algorithm works 
quite satisfactorily in this case also in that the computed V'4 is indeed constant 
as a function of the u-coordinate along J. Therefore, we show in Fig. 11 only 
a time profile for constant u. The line indicates the exact function while the 
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markers indicate the computed values. The relative error in this calculation 
(200 by 200 points) is a few percent in the region where the influence of the 
singularity is not too strong. 

Let us now discuss some of the issues related to the Bondi energy- momentum. 
This is an unexpectedly complicated issue which, in addition, depends on the 
global topology of the space-time under consideration. The standard definition 
used here is from [3]: 



where the integration is over a cut of JT". As it stands the formula is only valid 
under rather stringent simplifying assumptions. It is assumed that the surfaces 
O = const are null even away from J. This implies that J is non diverging and 
that the spin-coefficient r' vanishes. If these assumptions are not made, then 
the news J\f acquires additional compensating terms. 

The function W which appears in (5.11) is a function with conformal weight 
+1 on the cut satisfying the conformally invariant second order elliptic equation 



Here, the 3^ is the conformally invariant "eth" operator introduced in [3]. For a 
more standard form of this equation, we refer to [11]. The purpose of solving the 
equation (5.12) is to select out of the super-translation subgroup of the asymp- 
totic symmetry group (the BMS group) the normal subgroup of translations 
which is used to generate the energy-momentum. In the special case, where the 
metric on the cut has been scaled to be the standard unit sphere metric and 
where the cut sits within in such a way that the spin-coefficient r vanishes, 
then the equation (5.12) has four linearly independent solutions which can be 
taken as the first four spherical harmonics Yjm, I = 0, 1. Note, that r = can 
always be achieved in the neighbourhood of a single cut, because it only involves 
parallel transport of along the null generators. However, given a system of 
cuts and an adapted spin frame, this condition cannot, in general, be main- 
tained. Unfortunately, this is the case for the cuts appearing in the numerical 
treatment as intersections of J7 with the constant time hypersurfaces. A more 
thorough discussion of the general spherical case is left to another paper. 

Here we want to focus on our immediate interest, namely obtaining a formula 
for the Bondi mass on cuts with toroidal topology. In that case, the BMS group 
has a completely different structure. This is reflected in the fact that equation 
(5.12) on a torus has only a one dimensional solution space as opposed to the four 
dimensions in the spherical case. This means that the translation subgroup is a 
one-dimensional subgroup of the BMS group. Therefore, on toroidal cuts, there 
does not exist a four-vector of energy-momentum, but only a "Bondi scalar" , 
which we call the Bondi mass. 

In order to compare the evolution of that scalar with time in our special 
case we observe that for the initial data obtained from the exact solutions A3 
and Wl, the cuts are spanned by Killing vectors. This implies that on the cut 
all field variables are constant. Hence, we may take W = const as a solution of 
(5.12) and since W has to be a conformal density of weight -|-1 we take W = VA, 




(5.11) 



BlW = 0. 



(5.12) 
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with A being the area of the cut. Thus, we end up with the formula 



TUB = - 



VA<f {v2-^(AA-ap')} 



Jc 



(5.13) 



Of course, due to the constancy of the integrand on the cut we could have written 
this formula without the integral. However, we implement the formula with the 
integral because it averages over the numerical inaccuracies present from the 
interpolation process. In Fig. 12 is shown the (normalised) Bondi-mass for the 
A3 solution, which of course should remain constant. Similarly, Fig. 13 presents 
the Bondi mass for the Wl solution. Again, the solid line is the exact profile 
while the markers are the values obtained from the numerical solution. 

6 Conclusion 

In this article we have presented and discussed several issues concerning the 
numerical solution of the evolution part of the hyperboloidal initial value prob- 
lem for the vacuum conformal field equations. We have described a special 
case where the unphysical assumption was made that there exists a hypcrsur- 
face orthogonal Killing vector with closed orbits and no fixed points. This 
does not alter the essential issues. The numerical evolution scheme is a simple 
two-dimensional implementation of the well-known Lax-Wendroff method. The 
outer boundary is evolved using a stable eigen-field method. We have discussed 
various lapse choices and the features which appear when one specifies the gauge 
source function F as a function of the field variables. We have found a special 
choice for the shift vector which originates in the conformal properties of the 
system. This shift allows us to freeze null infinity on the grid while still leaving 
the usual freedom for specifying a shift vector in the interior. Finally, we have 
described how to obtain the local radiative information by simply "reading it 
off" J and transforming to the appropriate asymptotic spin frame. The global 
quantities like Bondi four-momentum are more difficult to determine and they 
are very different in our present case from the physical case where J has spher- 
ical sections. We have tested the code and the radiation extraction algorithm 
using exact solutions. We obtained good agreement between the analytical and 
the numerical solution. Unfortunately, the used exact solutions have an addi- 
tional Killing vector which makes them rather special even though we try to 
compensate for this by "warping" the coordinate system. Naturally, the next 
step has to be to obtain more general initial data. 

A The exact solutions 

We have used several exact solutions for numerical tests. Apart from the trivial 
ones which are simply Minkowski space in disguise, i.e., rescaled with an arbi- 
trary conformal factor there is the class of vacuum space-times with toroidal null 
infinities which have been constructed by Schmidt [2] for exactly that purpose. 
They are characterised by a solution w of a two-dimensional wave equation and 
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are defined as follows 

n = l{t'-z') (A.i) 

Given a solution of the two-dimensional wave equation [12] 

{t^ - z'^) {wtt - w,,) - 2t {3z^ + t^) wt - 2z {3f + z^) = 0, (A.3) 

one can obtain the function n{t, z) by quadratures. The coordinates x and y are 
Killing coordinates, each taking values in R. We identify the points {x, y) and 
(a;+ 1, y + T") to obtain the toroidal topology. In our applications, we always have 
r = 1. The simplest solutions of this type arc the ones obtained by choosing 
w{t, z) = (with n(t, z) — 0) and z) — A{t^ — v"^) for some constant a (with 
n{t,z) = -\A^{t^ + z'^Y). Note, that A = in the latter solution gives the 
former. The physical metric which corresponds to the first of these appears in 
the classification by Ehlers and Kundt [13] under the name A3 as the analogue 
of the Schwarzschild metric in plane symmetry. Here are the explicit expressions 
for the variables we use in the code. 



1 1 1 .2,,2 „ t (4A2[/2-3) 1 



2rr2 



V2VU ' '° V2 C/3/4 

^''~V^7U' ' ""'"'^ IPT' ' ' 

r _ * .^/TjAa'u' k * 4A^lP_UAU^ i^u2 

C,o-^VUe^ , ^4. = ^ e= , 

.V 1 1.2^,2 , 1 4A2[/2 - 12A;7+ 1 .2.,2 

^^° = *71c^'' ' "^" = 6 vu ' ' 

74r = -^^^^e^-^-^ = ^ ^^^^^ + 4 + 1 ^,2,2^ 

T, , t A^u' ^SU^A'^ + UA^U^ + GAU-S ^2^2 

Vu u^^^ 



1 ,2^2 /- %U^A^~AA^U^ + QAU + l .2^2 

4A2C/2 + 1 ,2,,2 



E = ^^e^ ^^^^ ct> = -2 '^T^^ e^^^^ 

2v2 y/U 

^ . V if— 1^2jj2 

^^°=-^^^^^ ' 

with U = t'^ + z^. All other functions either vanish or they are complex conju- 
gates of functions given above. 
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Figure 2: The local geometry for the timestep criterion 
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Figure 3: The proper time r (left) and the lapse N (right) in the "natural" 
gauge. The extremal values are Tj^jn = 1.831, Tmax = 2.988 and -/Vj^^jj^ = 0.315, 
-^max = 1.151. The black contours show the locations of the two J's. 




Figure 4: The proper time r (left) and the lapse N (right) in the harmonic gauge. 
The extremal values are t^^^ = 2.747, Tmax = 3.537 and -^V^nin = 0.03964, 
Nmax. = 0.1624. The zig-zag behaviour in the figure for A'' is due to a lack of 
sufficient resolution for the shading process. 
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Figure 6: Profiles of the lapse N for a fixed value of the Killing coordinate x 
and different values of p. Note the logarithmic scale. 
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Figure 7: The curvature invariant / along 2; = as a function of proper time. 
The line is the exact function, the points indicate computed values. 
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Figure 9: The proper time t (left) and the lapse N (right) in harmonic gauge 
with J7 freezing. The extremal values are Tj-^j^ = 1.831, Tmax = 2.988 and 
-^min ~ 0.315, A^max = 1.151. The black contours show the locations of the 
two J's. 
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Figure 12: The Bondi mass for the A3 solution normahsed against its initial 
value. 




Figure 13: The Bondi mass for the Wl solution normalised against its initial 
value. 



